Equilibrium insertion of nanoscale objects into phospholipid bilayers 
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Certain membrane proteins, peptides, nanoparticles and nanotubes have rigid structure and fixed 
shape. They are often viewed as spheres and cyhnders with certain surface properties. Single Chain 
Mean Field theory is used to model the equilibrium insertion of nanoscale spheres and rods into the 
phospholipid bilayer. The equilibrium structures and the resulting free energies of the nano-objects 
in the bilayer allow to distinguish different orientations in the bilayer and estimate the energy barrier 
of insertion. 



I. INTRODUCTION 

Biological membranes protect cells against the environ- 
ment and also provide for selective transport of nanoscale 
objects across the membranes^. Peptides, nanotubes, 
nanoparticles and other nano-objects with fixed shape 
can interact with phospholipid bilayers^, accumulate in- 
side the membranes'^ and penetrate into cells^"— . Surface 
activity of these objects raises issues of biocompatibil- 
ity and cytotoxicity^i^, which may limit their biomedical 
applications. Understanding and controlling the inter- 
actions of nanoscale objects with membranes is there- 
fore a key to the design of novel medical treatments and 
cell-active substances which can modulate the thermody- 
namic properties of cell membranes. 

Although membrane active peptides, membrane pro- 
teins, nanotubes and nanoparticles may significantly dif- 
fer in composition and chemical structure, the membrane 
activity implies the presence of certain common features 
providing for such activity. Phospholipid bilayers have 
alternating hydrophilic and hydrophobic regions which 
are characterized by the corresponding thickness. Thus, 
the size and the shape of nanoscale objects need to be 
compatible with the bilayer structure to ensure their in- 
sertion into the bilayer {e.g. the concept of hydrophobic 
matching^), cell penetrationiSiii or pore formatio n^^'^'^ . 
The thickness of the hydrophobic core defines the length 
scale for the size while the planar geometry of the bi- 
layer defines the template for the shape of nanoscale ob- 
jects able to interact or insert into phospholipid bilayers. 
The effect of the shape of nanoparticles interacting with 
phospholipid bilayers have been studied in Reffi^. The 
nanoparticles-lipid bilayer interaction is determined by 
the contact area and the local curvature of the nanoparti- 
cle at the contact point with the bilayer. Thus, nanoscale 
objects of different nature can be viewed in the first ap- 
proximation as geometrical objects described only by the 
shape and the surface propertieai^"— . 

Many theoretical studies of interaction of phospholipid 
bilayers with nanoparticles and even peptides consider 
the simplest shapes: a sphere^ii2ri2i and a cylinder— i22ri2S 
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FIG. 1. Three geometries of nanoscale objects considered in 
the text: insertion of a sphere (left), insertion of a perpen- 
dicular cylinder (middle), and insertion of a parallel cylinder 
(right) into the bilayer. 



with uniform or patterned surface properties. This allows 
to describe a broad range of objects of similar size within 
a single concept and provide a general picture on the 
mechanism and the molecular details of interaction with 
the bilayers. Coarse-grained Molecular Dynamics (MD) 
provides the most detailed information about the molecu- 
lar structure and the dynamics of interaction of nanoscale 
objects with phospholipid bilayer o'^d^ds, 19,23, 27- ^^ How- 
ever, this method becomes computationally expensive 
and time consuming when the size of the system and 
the number of interacting molecules increases^. Cal- 
culation of equilibrium properties and equilibrium free 
energies requires equilibration process which may exceed 
microseconds^. One of the alternatives is the use of 
the self-consistent field theories, where the equilibrium 
structure of the bilayer is found from the solution of the 
self-consistency equations. Recently, Single Chain Mean 
Field (SCMF) theory has been applied for modeling of 
the phospholipid bilayers at different levels of coarse- 
grainingS^. It was shown, that essential equilibrium prop- 
erties of the phospholipid bilayer such as the thickness, 
position of different groups, the compressibility constant 
and the area per lipid, can be reproduced with high accu- 
racy. High speed of calculations is obtained by the decou- 
pling of the correlations between the molecules and the 
fluctuations, taking into account the symmetries in the 
system. In this manner, the multi-molecule problem is re- 
duced to a single molecule in the external self-consistent 
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R = 4.05 A 




FIG. 2. The 3-beads model of the DMPC phospholipid 
molecule. R is the radius of the hydrophobic T-beads (red) 
and hydrophilic H-bead (blue), L is the bonds length, Rtt 
and Rhs are the interaction radii of the corresponding step- 
wise attractive potentials. 



field problem. The molecule conformations and the fields 
are found from the self-consistency condition and gives 
the free energy of the equilibrium structures as a result of 
calculations. This method has been applied to estimate 
the free energies of insertion of a carbon nanotube into 
a phospholipid bilayer^^ and to elucidate the patterning 
of a carbon nanotube for physical translocation through 
the bilayer— . 

In this publication, we first present equations of the 
SCMF theory for the interaction of phospholipid bilayer 
and nanoscale objects. This theory is then applied for 
calculation of equilibrium structures and free energies of 
resulting structures of phospholipids and nanoscale ob- 
jects, modeled as spheres and cylinders with given surface 
properties (Fig. [T|). 



II. THEORY 

The SCMF theory uses a coarse-grained representa- 
tion of a phospholipid molecule as a set of spherical 
beads interacting with the fields through simple effective 
potentials^. The simplest 3-beads model (Fig. ^ is able 
to reproduce with high accuracy the equilibrium struc- 
ture and properties of DMPC phospholipid bilayer, such 
as the average interfacial area per lipid, the thickness of 
the bilayer and the core, and the stretching modulus^^. 
In this model the phospholipid molecule is represented 
by two hydrophobic beads (T) and one hydrophilic bead 
(H) of the same radius 4.05 A, joined consequently by 
rigid bonds of length 10.0 A, and able to bend around 
the central T-bead. Two T-beads, belonging to differ- 
ent molecules, interact with energy gain ett = —2.10 kT 
at distances between their centers smaller than 12.15 A. 
Additionally, the H-beads interact with implicit solvent 
molecules, assumed as well to be represented by spher- 
ical beads of the same radius 4.05 A, with energy gain 
CHS = —0.15 kT at the distances closer than 12.15 A. 

The free energy of a system, containing N lipid 
molecules and the solvent, depends in general on the to- 
tal probability distribution function p of the system. In 
the mean field approximation one can split this multi- 



molecule probability distribution function into a product 
of the corresponding single-molecule distribution func- 
tions. This allows to write the total free energy of the 
3-beads model as a sum 



F = -SL-Ss + Un + Utt + Uhs + U, 



(1) 



where Sl and Ss are the entropy contributions of the 
lipids and the solvent, Uq describes the intra- molecular 
and position-dependent energies of the lipids, Utt and 
Uhs represent the inter-molecular interactions between 
the pairs of T-beads and between H-beads and the sol- 
vent, and the last term, Usteric, takes into account the 
hard-core steric repulsions between the molecules. 
The entropy part of the free energy is written as 

^Sl-Ss = N (In [pij)NA]) + J cs{r) In [cs{r)As] dr 

where ^(7) is the density of probability of a single lipid 
molecule to be in the conformation 7 (along with rela- 
tive positions of the beads in the molecule with respect 
to each other, the conformations also include their spa- 
tial position and orientation in the system), cs{r) is the 
solvent concentration at the point r, A and A5 are the de 
Broglie lengths for the lipids and the solvent respectively, 
which shift the absolute value of the free-energy but have 
no effect on the final results. The angular brackets de- 
note the average over ^(7), and the integral is taken over 
all the accessible volume of the system. 

The intra-molecular interaction energy Uq is simply the 
average over p(7) 



Uo = N{Uoh)) 



(3) 



where Uq (7) is the internal energy of the lipid molecule in 
the conformation 7. It can include interactions between 
the beads composing the molecules and interaction with 
the external fields and the walls (or impermeable immo- 
bile objects). 

The combined inter-molecular energy terms are writ- 
ten as 

N(N -1) f 

Utt + Uhs = / {uTh,r)) {cT{-f,r))dr{A) 



+N J {usij,r)) cs{r)dr 

where UT{j,r) and us{'j,r) are the interaction potentials 
for T-beads and the solvent at the point r, created by the 
lipid molecule in the conformation 7, and ct(7,^') is the 
concentration of T-beads of the lipid in conformation 7 
at the point r. 

[5] is valid only if the potentials ut(7,?") and 1*3(7, r) 
are finite and smooth functions of r and therefore cannot 
contain any hard-core steric repulsion interactions, which 
are included in a separate term, Usteric, representing the 
incompressibility condition at every point r of the system 



U 



steric 



X{r)i(j)o-vsCsir)-N{^ij,r)))dr (5) 
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Here A(r) is the Lagrange multiplier and (f>Q is the to- 
tal volume fraction occupied by the lipids and the sol- 
vent {0 < (f)Q < 1), vs denotes the volume of the solvent 
molecule, and (f>{j, r) is the volume fraction at the point 
r occupied by the lipid in the conformation 7. 

To find the equilibrium distribution, one should mini- 
mize the total free energy ([IJ with respect to p{'^) and the 
solvent distribution csir) subject to the incompressibil- 
ity condition ([S]) . The minimization leads to the following 
set of integral equations 



Similarly, the integrals over the space are replaced by the 
sums as 



P(0 = exp 



-U„{0-{N-1) / UT{^,r){cTh,r))dr- 



us{^,r)cs{r)dr + / A(r)(/)(^,r)dr 



■ysA(r) = In [vscs{r)] + N (1*5(7, r)) 



vscs{r) = (l)s{r) = 0o -^(0(7,^)) 



(6) 



(7) 



(8) 



where the normalization factor Z ensures that 
/ p{'-f)dj = 1. These equations can be solved nu- 
merically, if the integrals and the averages are replaced 
by appropriate sums, which reduce the problem to 
solution of algebraic equations. To discretize the space, 
the simulation box is divided into a set of auxiliary cells, 
where all spatial dependent variables have a constant 
value. The division into cells should also take into 
account the symmetry of the solution, which provides 
the speed up of the calculations. For example, if the lipid 
system is expected to self-organize into a flat bilayer, the 
division of the simulation box according to ID geometry 
into parallel slits results in considerably smaller number 
of cells than the 3D division into a cubic grid with the 
same resolution. The conformational space 7 of a single 
lipid is sampled with the Rosenbluth algorithm^^, the 
conformations are placed in the box and the spatial 
distributions Mr (7,7'), Ms (7,7'), ct(7,'"), 4>{^,r) as well 
as other necessary conformational-dependent properties 
are calculated. Finally, the averages of the fields / in 
the equations [B]|H] are replaced by the sums 



(/(7)) = -E 



J fir)dr = J2m 



(10) 



(9) 



where fi is the value of the field in the i-th auxiliary cell 
of the simulation box and Vi is the volume of the i-th 
cell. As a result, the system of integral equations [niS] is 
reduced to the closed set of algebraic equations of finite 
but huge number n of unknown variables p-y. Since the 
larger k, the more accurate the solution is, the practical 
choice of n which provide a reasonable accuracy is lim- 
ited by the representative sampling. The easiest way to 
check the accuracy is to perform several simulations with 
independently generated random samplings and compare 
their results. Thus we start our simulations with some 
small value of k, and then perform series of similar simu- 
lations for larger k till the moment when further increase 
of K will not lead to significant change of calculated re- 
sults. For the simulations reported in this work we used 
K = 750000, which provide accuracy of the calculated 
energies about ±10 kT. 

The number of auxiliary cells in the box should be 
high enough to provide acceptable spatial resolution, 
but if it is relatively small one can decrease the num- 
ber of independent variables by introducing new in- 
dependent variables, average concentration of T-beads, 
CT.i = {cT{'y,ri)), potential of the solvent, us,i = 
{usi^jVi)), and the volume fraction occupied by lipids, 

= {4>{lJ' i)) ■ Depending on the method of solution of 
equationtlj, such reduction of the number of independent 
variables can speed up calculations. 

The equations of the SCMF theory are written for the 
canonical (NVT) ensemble. However, the simulation of 
the phospholipid bilayer with inserted nanoscale object 
would require the grand-canonical ensemble, since the 
inserted object may expel the phospholipids out of the 
simulation box. Fortunately, the SCMF technique is fast 
enough to perform series of canonical calculations for a 
part of bilayer with variable number of lipids N in the 
simulation box. One can consider then the simulation 
box as an open part of a larger closed system and thus 
estimate the free energy F* of the whole system as a 
function of N . The equilibrium value of N will corre- 
spond to the minimum of -F*, thus series of canonical 
SCMF calculations are equivalent to a grand-canonical 



where is the value of this field, and is the proba- 
bility of the conformation 7 belonging to the generated 
sampling (thus p-^ is determined only over the sampling, 
while ^(7) used above was determined in the full confor- 
mational space of the molecule), is the value of the 
Rosenbluth weight [i.e. biased probability of the con- 
formation 7 is introduced into the sampling during the 
generation), k, is the size of the conformational sampling. 



^ According to Ref,^^, the independent variables are cy ^, cjj i 
and 4>i. One can show that the expression Af (115(7, r)), which 
is the average potential experienced by the solvent molecule at 
the point r, is approximately equal to Nes (c//(7,r)}, where ts 
is the integral of the H-S interaction potential around a solvent 
molecule. One can easily replace the variables itg ^ by cjj This 
was indeed necessary in the general case, considered in the pre- 
vious article, but in this context it can be avoided for the sake 
of simplicity. 
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calculation. In the simulations performed in the present 
work (size of the system about 100.0 x 100.0 x 60.0 A, 
two-dimensional geometry with spatial resolution 2.0 A, 
and sampling size k = 750000), calculation of one energy 
point in the grand-canonical ensemble with current ver- 
sion of our SCMF code, demanded, on average, about 
one day on a 32-cores machine. Full three-dimensional 
simulations should be significantly slower, but we believe 
that further improvements of our code will significantly 
speed-up its' performance, thus making possible three- 
dimensional simulations in a reasonable time. 



III. RESULTS AND DISCUSSION 

The equations of the previous section allow to calcu- 
late the equilibrium properties and the free energy of the 
phospholipid bilayer and its interaction with nanoscale 
objects with a fixed shape. Simulation of a phospholipid 
bilayer is a one dimensional problem if we neglect the 
bilayer bending. Since the free energy of a homogenous 
bilayer is proportional to the area (or the number of lipids 
in the bilayer) and the interfacial area per lipid and the 
bilayer thickness are constant along the bilayer, the free 
energy of a large system F* can be written as 



p* 



V*fs+N*f{A), f{A) 



F{V,N,A)-Vfs 
N 



fs 
(11) 

where V* is the volume and N* is the total number of 
lipids of the large system, fs is the free energy density of 
the region occupied by the pure solvent, and f{A) is the 
free energy per lipid of the bilayer, which depends on the 
average surface area per lipid A. 

The result of the SCMF simulation is the free energy 
F{V, N, A) of the box of the total volume V with the 
number of lipids N. These expressions allow to con- 
struct the free energy per lipid f{A). The equilibrium 
area per lipid ^o, which minimizes the total free energy 
of a large system F* , can be found by the minimization 
of f{A), obtained from the series of SCMF simulations, 
while the simulation corresponding to Aq gives the equi- 
librium distribution of T- and H-beads inside the bilayer 
and the equilibrium thickness (Fig. |3]). 

Similar approach is applied to study the insertion of a 
rigid spherical or cylindrical nano-object into the phos- 
pholipid bilayer. A part of the object located inside 
the simulation box is represented by a region prohibited 
for lipids and the solvent to enter, and the surface of 
the object can interact preferentially with phospholipid 
beads. This restricted area is taken into account during 
the sampling generation by the additional bias and the 
corresponding contribution to the intra-molecular energy 
Uo^j). In the case of a spherical object and a cylindri- 
cal object oriented perpendicular to the bilayer one can 
use the two-dimensional cylindrical geometry with the 
symmetry axis passing through the center of the sphere 
(cylinder) perpendicular to the bilayer plane. The sim- 
ulation of a cylindrical object in the parallel orientation 
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FIG. 3. Mean field snapshot of DMPC bilayer at equilibrium 
within the 3-beads model (left), and the corresponding con- 
centration profiles (right) of hydrophobic (red), hydrophilic 
(blue) and total (black) volume fractions of the beads. 



has been done within two-dimensional grid geometry i.e. 
assuming infinite cylinder. To limit the bending of the 
bilayer, which permits the inserted object to escape from 
the bilayer, the hard walls on the top and bottom of the 

(^niu]0(j)ion box have been introduced. 

yg As%i^ming that a part of the membrane outside the 
simulation box is not perturbed by the insertion of the 
object and, thus, is in equilibrium state, we estimate the 
total free energy change due to the insertion of the object 
to a given position p in the box as 



AF* = F{V, N,p) - Nf{Ao) - {V - V{p))fs (12) 

where F{V, N,p) is the free energy of the box and V{p) 
is the volume occupied by the object inside the box. 

The results of the free energy calculations for spheres 
and cylinders of diameter 2.43 nm with homogeneous sur- 
face properties are summarized in Fig. |4l The surface 
properties are described by a single interaction parameter 
between the T-beads and the surface, e. One can clearly 
see, that the free energy changes due to the insertion of a 
sphere, or a cylinder in perpendicular orientation which 
has a similar contact area with lipids, have relatively close 
values for all considered levels of hydrophobicity e. Thus, 
the shape of the object is not so important as the con- 
tact area with lipids. This conclusion is consistent with 
the previous observational^. In contrast, the cylinder in 
the parallel orientation has a larger contact area with 
phospholipid tails and disturbs much larger part of the 
bilayer. Thus, the parallel orientation has sufficiently 
lower equilibrium free energy than the perpendicular ori- 
entation (the points in the plot correspond to a cylinder 
with the length 10.0 nm, while the numbers in the insets 
correspond to the energy per nm). 

Hydrophilic and hydrophobic objects induce different 
rearrangement of lipids in the bilayer (see the insets of 



cylinder 18.2kT/nm cylinder -1.9kT/nm 



sphere -301 .5 kT 




FIG. 4. Snapshots of the equihbrium insertion of a spherical and a cylindrical object into the DMPC bilayer. The color of 
the objects in the insets illustrates the hydrophobicity level (interaction with T-beads): e — 0.0 (blue), e = —6.3 kT (purple), 
e = —10.0 kT (red). Parallel orientation of a cylinder (denoted by the energy/nm), corresponds to the cylinder of length 10.0 



Fig. 13]). Hydrophilic objects of similar size create pores 
of similar size with the heads of the lipids oriented inside 
the pores. Hydrophobic objects attract the tails of phos- 
pholipids and they tend to come closer to the surface. 
For most hydrophobic objects, the meniscus provoked by 
the wetting of the hydrophobic object can be formed. 
For intermediate interaction parameters, the pore is not 
complete and is combined with the partial wetting at 
the edges, thus provoking the thinning of the bilayer in 
contact with the nano-object. This diagram also sug- 
gest the preferential orientation of cylindrical objects in 
the bilayer. Hydrophilic cylinder has lower free energy 
in perpendicular orientation, while the most hydropho- 
bic cylinder has lower free energy in parallel orientation. 
For an intermediate hydrophobicity, there is a transition 
in orientation from perpendicular to parallel orientation 
which depends also on the length of the cylinder. For 
e = —6.3 kT the transition from perpendicular to paral- 
lel orientation occurs at the length of the cylinder 10 nm, 
where both orientations have almost the same energies. 
Longer cylinders would favor parallel orientation. 

The uptake of a hydrophobic particle into the core of 
a bilayer is not the only equilibrium solution of the equa- 
tions. There exists a solution which corresponds to a 



hydrophobic particle covered by phospholipids floating 
around and coexisting with a bilayer (third snapshot from 
the top in the right column of Fig. 0]). We have compared 
the free energy of such structure with the energy of the 
inserted sphere into the bilayer in the most hydrophobic 
case, e = —10.0 kT. The bilayer with the incorporated 
sphere has a lower energy, while the free energy difference 
with the sphere in the solution is of the order 50.0 kT, 
indicating that the hydrophobic sphere would unlikely 
escape. Similar conclusion can be drawn for the case of 
cylindrical particles lying parallel to the bilayer. 



IV. CONCLUSIONS 

We have developed the SCMF theory for the inter- 
action of nanoscale objects with phospholipid bilayers. 
This numerical method provides for detailed microscopic 
information about the structure of phospholipid bilayers 
and their essential equilibrium thermodynamic proper- 
ties. The method gives the total equilibrium free en- 
ergy of self-assembled structures as the output of calcu- 
lations, which allows for direct comparison of different 
equilibrium structures involving nanoscale objects and 
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phospholipid bilaycrs. 

Wc have apphed the SCMF theory to investigate 
the interaction of phosphohpid bilayer with cylindrical 
and spherical nanoparticles of diameter 2.43 nm and 
surface properties ranging from hydrophilic to strongly 
hydrophobic. We have shown that the shape of the 
nanoparticles is not as important as the contact area with 
the lipids. Hydrophobic and hydrophilic objects induce 
different reorganization of lipids around the nanoparti- 
cles. More complex nano-objects will be considered in 



future works. 
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